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Abstract 

The natural convection of incompressible fluid in a porous medium causes for some 
boundary conditions a strong non-uniqueness in the form of a continuous family of 
steady states. We are interested in the situation when these boundary conditions 
are violated. The resulting destruction of the family of steady states is studied via 
computer experiments based on a mimetic finite-difference approach. Convection 
in a rectangular enclosure is considered under different perturbations of boundary 
conditions (heat sources, infiltration). Two scenario of the family of equilibria are 
found: the transformation to a limit cycle and the formation of isolated convective 
patterns. 
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1 Introduction 



The Darcy model of two-dimensional incompressible fluid convection in a 
porous medium [13,3,15,16] demonstrates the nontrivial effect of the appear- 
ance of a continuous family of steady states (equilibria) with a variable spec- 
trum of stability. It was explained in [15] that the cosymmetry can lead to this 
phenomenon. Cosymmetry differs from symmetry due to variability of spec- 
trum where all members of the family have the identical spectrum of stability. 
The effects of cosymmetry in Darcy convection was analyzed theoretically in 
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[15,16,17]. The computation of the families of steady states and their evolution 
was carried out in [4,7,5,6,8,10]. Partial instability of the family (the occur- 
rence of unstable steady states) was studied in [4,5,6]. The effects of collision 
and reorganization of the families were analyzed in [5,6,8]. 

A related interesting phenomenon is the destruction of a family of equilibria 
under some perturbations [2,17,12]. In [2] weak seeping of the fluid through the 
boundaries was considered. The breakdown of the family of stationary solu- 
tions was analyzed using the normal forms technique. V. Yudovich [17] studied 
the violation of the cosymmetry under low-intense heat sources or weak fluid 
infiltration and found that the continuous family of equilibria disintegrates or 
even disappears. It was predicted that in the case of heat sources a family of 
steady states may disintegrate to produce several separate stationary regimes. 
The case of weak infiltration may lead to the appearance of a long-period limit 
cycle instead of the family of equilibria [2,17]. 

In this paper we study the problem of the destruction of the family of steady 
states in the planar problem of Darcy convection in a rectangle. The finite- 
difference approximation of the Darcy equation allows us to consider the per- 
turbation listed in [17] but our numerical studies are not restricted to the 
cases of small disturbances of heat source and low velocity of infiltration. We 
carried out several computer experiments with the harmonic perturbations of 
the temperature on the bottom (the case of a finite number of heat sources 
located on the boundary). The disintegration of the family of continuous equi- 
libria is demonstrated numerically by showing the appearance of limit cycles 
and the formation of isolated convective patterns. 

The paper is organized as follows. In Section 2 we give the Darcy equation and 
the cosymmetry property briefly. In Section 3 the finite-difference discretiza- 
tion is outlined. In Section 4 we present the numerical results for different 
perturbations of the boundary conditions. 



2 Darcy model and cosymmetry 

We consider an incompressible fluid which saturates a porous medium in a 
rectangular enclosure V = [0, a] x [0, b] which is heated below. The dimension- 
less equations of the fluid flow in a porous medium according to Darcy law 
can be written as [16]: 



6 t = A9 + \ij x + J(9,4>) =F 
= Aip - 9 X = G, 



(1) 
(2) 
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where A is the Laplacian and J(8, tji) denotes the Jacobian operator over 
(x,y) J(8,ip) = 6 x ip y — 6yip x . The unknown functions ip(x,y,t) and 9(x,y,t) 
denote, respectively, the stream function and the perturbation of temperature 
from the basic state of rest with a linear conductive profile, A is the Rayleigh 
number. 

The boundary conditions are: 



tp(x,y,t) = fiy, (x,y)edV, 

here the function f(x) gives the nonhomogeneous distribution of the temper- 
ature at the bottom with magnitude 77, and fi is the velocity of infiltration 
through the lateral sides of an enclosure. 

The problem is supplied by the initial conditions 



where 6q denotes the initial perturbation of the temperature from the linear 
profile. For a given 9q, the stream function can be obtained from (l)-(3) as 
the solution of the Dirichlet problem via Green's operator. 

When 77 = fi = the system (l)-(3) imposes the equilibrium 8 = ip = 0. 
This state of rest remains stable while A < An. The critical values \j are 
given by \j = n 2 {i 2 /a 2 + j 2 /b 2 ), i,j G Z. It was shown in [16] that the first 
critical value An is two-times for any domain V and larger values of A implies 
an onset of the family of steady states. This result is a consequence of the 
cosymmetry that exists for the system (l)-(3) when 77 = fi = [15]. Indeed, 
the vector-function (8, —ip) is orthogonal to the right-hand side of (1), (2) in 
L 2 , i.e Jx,(Fip — G8)dxdy = 0. In the case of nonzero 77 and ji we obtain using 
integration by parts and Green's formulae 



0(x, 0, t) = rjf{x); 8(x, y, t) = 0, (x, y) G &D \{y = 0}; 



(3) 



8(x,y,0) = 9 (x,y), 



(4) 



a 



a 




77 / ip y (x, 0, t)f(x)dx — fib / 9 y (x,b,t)dx 



v 











b 








Thus, the cosymmetry condition fails whenever fi 7^ or 77 7^ 0. 
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3 Finite difference discretization 



We use here the discretization and numerical scheme derived in [9]. Let us 
introduce the uniform arrangement of the nodes on the interval for spatial 
coordinates: 



x n = nh x , n — 0, • • • , N + 1, h x = a/(N + 1), 
y m = mh y , m = 0,---,M + l, h y = b/(M + l). 

The temperature 9 and stream function ip are defined at the nodes (x n , y m ). It 
is useful to introduce the auxiliary (staggered) nodes on each of the coordinates 
Xn+i/2 = \{x n + x n+1 ), n = 0,---,N, and y m+ i/ 2 = \{y m + Um+i), m = 
0,-..,M. 

To approximate equations (l)-(2) we introduce difference and averaging op- 
erators based on the two-node stencils 



/ r r\ _ fn+l,m fn,m / <• «\ _ fn,m+l fn,m /r\ 

\0xJ )n+l/2,m ~ ^ , \0y J )n,m+l/2 ~ ^ , {3) 

/ s f\ _ fn+l,m fn,rn ( X f\ _ fn,m+l + fn ,rn / n \ 

\OQxJ )n+l/2,m ~ ~ , \°0yJ )n,m+l/2 ~ ~ • {") 

These formulas are valid for integer and half-integer values of n and m. Using 
operators (5)-(6) we derive the operators on three-node stencils 

D x = S x 5 0x = S 0x 5 x , D y f = S y S 0y = S 0y 5 y , (7) 



the averaging and difference operators on a four-point stencil 

do = 5o x 5o y = 5o y 5o x , d x = 5 x 5o y = 5o y S x , d y = S y So x = So x S y , (8) 

and the discrete analog of Laplacian 

(A,/)„, m = (5JJ + 5 y 5yf) nm . (9) 

To approximate the Jacobian we use the Arakawa scheme [1] 



Jn,m(W)= \[D x {9Dy^)-D y {9D x ^)} nm (10) 
2 

+ g \d x (d 9d y i)) - d y (d 9d x i)) ] n m . 
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The discretized boundary conditions are formulated as follows 

^n,o = 0, e nt0 = rjf{x n ), n= 0,...,N + 1, (11) 
VVi,m+i = 0, 9 n , M +i = 0, n= 0, ...,iV + l, 

#0,™ = N +l,m = 0, tl>0, m = lj>N+l,m = Wm, ™ = 0, . . . , M + 1. 

Using operators (7)-(10) we obtain the semi-discretized form of (1) and (2): 

= ^nm = (A h + AZ^nm + (12) 

= G nm = (A h il>-D x 9) nm , (13) 

here n = 1, • • • , N, m = 1, • • • , M, the dot denotes a differentiation on t. 

After discretization the resulting system of ordinary differential equations 
(12)-(13) may be written in the following form 

= (A + \BA~ 1 B)Q + F{Q) = $(9, A, n, rj). (14) 

Here 6 is vector of unknowns (temperature), matrices A and B are the ap- 
proximation, respectively, of the Laplacian A and the first derivative on x, 
where F(Q) represents the nonlinear Jacobian term. 

To find nonstationary regimes or to converge to stable equilibrium we apply 
the Runge-Kutta method of fourth order. To compute a family of equilibria 
in the system (14), we apply the technique [7] based on the cosymmetric 
version of the implicit function theorem [16] . Different realization for ordinary 
differential equations was also developed in [4]. 

When A is slightly larger than An, then all points of the family are stable 
[16]. Starting from the vicinity of unstable zero equilibrium we integrate the 
ordinary differential equations (14) by the Runge-Kutta method up to a point 
6 close to some stable equilibrium on the family. Then we correct the point 6 
by the Newton method. To predict the next point on the family we determine 
the kernel of the linearization matrix (Jacobi matrix) at the point G>o and use 
the Adams-Bashford method. We repeat this procedure to obtain the entire 
family. 



4 Numerical results 

Here we mainly consider convection for a wide enclosure with a = 2 and 6=1. 
Our previous computation of the families of steady states [7,8] showed that a 
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sufficiently good description of a number of the transitions can be obtained 
on a rather coarse grid. Several hundreds of internal nodes was enough to 
describe an appearance of a primary family of steady states [7] and follow a 
number transformations [8]. So we used the discretization with 16 x 8 nodes 
and compared with results on 32 x 16 nodes to study the destruction of a 
family of steady states under inhomogeneity of boundary conditions. 

To present results we apply the special Nusselt coordinates [4] 

b a 

Nu h = J 9 x (^,y)dy, Nu v = J 9 y (x,0)dx, (15) 





where Nuh corresponds to the integral value of the heat flux from left to right 
defined for the centered vertical section. The value Nu v may be considered as 
a cumulative heat flux through the bottom of the container. 



4-1 Infiltration through the lateral sides 



In [2,17] a scenario of the transformation of a family of steady states to a 
limit cycle of the large period was described. Fig. 1 illustrates this scenario 




Fig. 1. Family of steady states (dashed) and limit cycles (solid) for several infiltra- 
tion velocities u: A = 60 (left) A = 240 (right) 

for two Rayleigh numbers in the case of small and large infiltration through 
the lateral sides. The case with A = 60 corresponds to the initial stage of 
the evolution of the family of equilibria after its birth. The Rayleigh number 
A = 240 shows the rather developed convective patterns. If infiltration velocity 
/i is small then the trajectory of a limit cycle lies in vicinity of the family of 
equilibria. When the velocity /i increases a limit cycle goes far away from the 
family. Moreover, large infiltration causes degradation of the nonstationary 
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convective regime. All fluid particles move horizontally without the formation 
of a convective pattern. With the increase of the Rayleigh parameter A the 
convective regime becomes able to prevent a pure infiltration. Computation 
at A = 240 shows that for small velocity \i a limit cycle occurs near a curve 
of the family of equilibria. Increasing of \x shifts the limit cycle to the range 
of small Nusselt values Nuu- Afterwards we see a deformation of a limit cycle 
and rather large velocity of infiltration leads to suppression of the convection. 

The period of a limit cycle T becomes smaller with increasing velocity fi 
and depends also on Rayleigh parameter A. One can see in Fig. 2 that the 
dependence T(/z) is practically the same for different values of A. 




Fig. 2. Dependence of the period of a limit cycle from the velocity of infiltration \x 
for several A 



4-2 Inhomogeneous heating from below 



For inhomogeneous heating and without infiltration (/i = 0) we found two 
scenarios of destruction of the family of equilibria. We consider the case of the 
harmonic distribution of the bottom temperature f(x) = sin^nnx / a) . It allows 
us to consider the case with zero average deviation from uniform heating. 
Through the experiment we change the intensity of heating (Rayleigh number 
A), the magnitude of inhomogeneity r\ and the amount of half-harmonics k. 
Depending on the parameters, the family of steady states may transform to a 
limit cycle or to a number of isolated equilibria. 

Computation with two harmonics in temperature distribution (k, = 4) shows 
that the family of steady states disintegrates onto two isolated equilibria, see 
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Fig. 3. Family of steady states and convergence from different initial points (circles) 
to the isolated stationary states (stars) under inhomogeneous heating below: r] = 4 
(dashed), r\ = —4 (dashed-point); A = 70, k = 4 

Fig. 3. One of these equilibria is stable while the second is unstable, it depends 
on the sign of magnitude r\. 

The computer experiment was started from points taken from the curve of 
the family of equilibria. For small values of magnitude 77 the trajectories were 
very close to the curve of the family of equilibria. In this case we found very 
slow motion to the isolated stationary state. 

In the case of one harmonics (k, = 2) the scenario of destruction of a family 
onto isolated equilibria was observed when A < 180. The transformation of 
the family of steady states to a limit cycle was found for A > 180. Fig. 4 gives 
relative location of the family of steady states and two limit cycles. For small 
magnitudes 77 the trajectory of the limit cycle lies close to the family. 

The convective patterns and their number depend on the number of half- 
harmonics k. Fig. 5 shows the results of disintegration of the families for two 
Rayleigh numbers: A = 120 and A = 180. We changed k and fixed a magni- 
tude of heating perturbation 77 = 8. For k = 2 and k = 4 only one isolated 
equilibrium was observed. The 6 produces two isolated stable equi- 

libria instead of the family of steady states. The unstable equilibria also exist 
and may be obtained using computation with 77 = —8. Here we see a change 
of stability properties for isolated equilibria to which the family of steady 
states transforms. Stream functions for realized convective patterns are given 
in Fig. 6. It is clearly seen that one state for k = 6 looks as if the equilibrium 
being formed under temperature distribution with two harmonics (k, = 4). 
The second possible state for k — 6 is rather different. These experiments 
show that we can control the desirable state by the appropriate choosing of 
magnitude and amount of harmonics. 
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Fig. 4. The limit cycle (dashed) and the family of steady states (solid); A = 190, 
rj = 8 (curve 1), rj = 16 (curve 2), k = 2 



Nu 



60 



-60 



-120 



^20 




Nu 



120 



-120 



-240 



330 



350 



370 



Nu, 




Fig. 5. Families of steady states (solid) and isolated states (stars) for heating with 
magnitude of perturbation rj = 8 and different number of half-harmonics k. Left: 
A = 120, k = 2 (1), k = 4 (2), k = 6 (3), k = 6 (4); Right: A = 180; k = 2 (5), k = 4 
(6), k = 6 (7), k = 6 (8) 

5 Conclusion 



We have studied convection in the rectangular enclosure under the following 
perturbations of boundary conditions: heat sources at the bottom and infil- 
tration through the lateral sides. Nonzero boundary values destroy the family 
of steady states which existed in the case of Dirichlet boundary conditions. 
Computer experiment allows to consider the case of large magnitudes of per- 
turbations. Two scenario were studied: transformation of the family to a limit 
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Fig. 6. Stream function of isolated steady states for different k: A = 120 (left), 
A = 180 (right); n = 8 

cycle and destruction of a family on a number of isolated convective patterns. 

We found through computation that convergence to the isolated steady states 
may be very slow. It looks as a memory effect when the system remembers the 
location of the former family of equilibria. It should be noted that the dynamics 
on finite time span is very close to the results on selection of equilibria under 
bottom heat sources [14]. This may be used to control the preferable convective 
state via variation of boundary conditions. 
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